Research on optical imaging techniques is ongoing in medical fields to acquire information on the inside of an object by irradiating light onto the object, and PAT (Photoacoustic Tomography) is one of these techniques. PAT utilizes the photoacoustic effect when a light absorber (a region having a high absorption coefficient) in the object absorbs energy of the irradiated light, expands its volume, generates a photoacoustic wave, and visualizes information related to the optical characteristic values inside the object. To visualize information, an acoustic wave is detected at a plurality of locations surrounding the object, and acquired signals are mathematically analyzed.
If a photoacoustic diagnostic apparatus based on the PAT technique is used, such information as initial sound pressure distribution and light energy absorption density distribution inside the object can be acquired, and the acquired information can be used for specifying a location of a malignant tumor that involves the growth of new blood vessels, for example. In the following, the description on initial sound pressure distribution includes description on light energy absorption density. Generating and displaying a three-dimensional reconstructed image based on such initial sound pressure distribution is useful in knowing the internal state of biological tissue for diagnosis.
On the other hand, along with the recent advancements of information processors and the increase in data capacity, three-dimensional images of the human body, such as CT and MRI, are used at higher frequencies in the medical field. For medical image diagnosis, it is desirable to save all the data for a long period of time in order to compare the test results of a plurality of modalities, and observe progress after surgery. However three-dimensional image data normally has large capacity, and redundant data must be minimized for long term data storage.
For the format of three-dimensional image data, a standard format, such as volume data format, is better than an application-specific format in terms of versatility which allows the use of various software, and easy data analysis. Therefore in a photoacoustic diagnostic apparatus as well, it is preferable to output a three-dimensional image in volume data format.
Now the characteristics of a reconstructed image by PAT will be described. In the case of PAT, if the time variation of an acoustic wave is measured at various points on a closed spatial surface (a spherical measurement surface in particular) that surrounds the entire object, using an ideal acoustic detector (wideband, point detection), the initial sound pressure distribution generated by photo-irradiation can theoretically be completely visualized. It is also known that the initial sound pressure distribution generated by photo-irradiation can be reproduced almost perfectly if columnar or planar measurement is performed on the object, even if a closed space is not used (see NPL 1).
Expression (1) is a partial differential equation called a “photoacoustic wave equation”, and by solving this equation, acoustic wave propagation from the initial sound pressure distribution can be described, and where and how the acoustic wave could be detected can be theoretically performed.
                    [                  Math          .                                          ⁢          1                ]                                                                                  (                                                            ∇                  2                                ⁢                                  -                                      1                                          c                      2                                                                                  ⁢                                                ∂                  2                                                  ∂                                      t                    2                                                                        )                    ⁢                      p            ⁡                          (                              r                ,                t                            )                                      =                              -                                          p                0                            ⁡                              (                r                )                                              ⁢                                    ∂                              δ                ⁡                                  (                  t                  )                                                                    ∂              t                                                          (        1        )            where r denotes a position, t denotes time, p(r, t) denotes time variation of the sound pressure, p0(r) denotes initial sound pressure distribution, and c denotes sound velocity. δ(t) denotes a delta function that represents the shape of the light pulse.
Reconstructing an image by PAT means deriving the initial sound pressure distribution p0(r) from the sound pressure pd (rd,t) acquired at a detection point, which in mathematics is called an “inverse problem”. The UBP (Universal Back Projection) method, which is a representative image reconstruction method based on PAT, will now be described. The inverse problem to determine p0(r) can be accurately solved by analyzing the photoacoustic wave equation of Expression (1) in the frequency space. UBP is this result expressed in the time space. Finally Expression (2) shown below is derived.
                    [                  Math          .                                          ⁢          2                ]                                                                                  p            0                    ⁡                      (            r            )                          =                              -                          2                              Ω                0                                              ⁢                      ∇                          ·                                                ∫                                      S                    0                                                                                                          ⁢                                                                            n                      ⋒                                        0                    S                                    ⁢                                                                          ⁢                                                                                    dS                        0                                            ⁡                                              [                                                                                                            p                              0                                                        ⁡                                                          (                                                                                                r                                  0                                                                ,                                t                                                            )                                                                                t                                                ]                                                                                    t                      =                                                                                                r                          -                                                      r                            0                                                                                                                                                                                                                              (        2        )            where Ω0 denotes a solid angle of the entire measurement area S0 at an arbitrary voxel (unit region).
This expression can be simplified and transformed into Expression (3) shown below.
                    [                  Math          .                                          ⁢          3                ]                                                                                  p            0                    ⁢                                          ⁢                      (            r            )                          =                              ∫                          Ω              0                                                                      ⁢                                    b              ⁡                              (                                                      r                    0                                    ,                                      t                    =                                                                                        r                        -                                                  r                          0                                                                                                                                          )                                      ⁢                                                  ⁢                                          d                ⁢                                                                  ⁢                                  Ω                  0                                                            Ω                0                                                                        (        3        )            where b(r0,t) denotes projection data, and dΩ0 denotes a solid angle of a detector dS0 to an arbitrary observation point P. The initial sound pressure distribution p0(r) can be acquired by performing back projection of this projection data according to the integration of Expression (3).
b(r0,t) and dΩ0 are given by Expression (4) and Expression (5) shown below.
                    [                  Math          .                                          ⁢          4                ]                                                                      b          ⁡                      (                                          r                0                            ,              t                        )                          =                              2            ⁢                                                  ⁢                          p              ⁡                              (                                                      r                    0                                    ,                  t                                )                                              -                      2            ⁢                                                  ⁢            t            ⁢                                          ∂                                  p                  ⁡                                      (                                                                  r                        0                                            ,                      t                                        )                                                                              ∂                t                                                                        (        4        )                                          ⅆ                      Ω            0                          =                                            dS              0                                                                    r                -                                  r                  0                                                                            ⁢          cos          ⁢                                          ⁢          θ                                    (        5        )            where θ is an angle formed by the detector and an arbitrary observation point P.
If the distance between the sound source and the measurement position is sufficiently long with respect to the level of the sound source (long distance sound field approximation), Expression (6) shown below is used.
                    [                  Math          .                                          ⁢          5                ]                                                                      p          ⁡                      (                                          r                0                            ,              t                        )                          ⪡                  t          ⁢                                    ∂                              p                ⁡                                  (                                                            r                      0                                        ,                    t                                    )                                                                    ∂              t                                                          (        6        )            
In this case, b(r0,t) is given by Expression (7) shown below.
                    [                  Math          .                                          ⁢          6                ]                                                                      b          ⁡                      (                                          r                0                            ,              t                        )                          =                              -            2                    ⁢          t          ⁢                                    ∂                              p                ⁡                                  (                                                            r                      0                                        ,                    t                                    )                                                                    ∂              t                                                          (        7        )            
Thus in the image reconstruction based on PAT, the projection data b(r0,t) is determined by time-differentiating the detection signal p(r0,t) acquired by the detector, and is back-projected according to Expression (3), whereby the initial sound pressure distribution p0(r) is determined (see NPL 1).
Expression (1), used for determining Expression (3), however, assumes “constant sound velocity”, “measurement from every direction”, “impulse type photo-excitation”, “acoustic wave detection in broadband”, “acoustic wave detection at a point” and “continuous acoustic wave sampling”. In reality it is not easy to implement an apparatus that satisfies these assumptions.
For example, it is actually difficult to acquire acoustic wave detection information on the total closed spatial surface surrounding the entire object. Furthermore, in order to increase the acoustic wave measurement region, the size and number of elements of the acoustic detector, and a signal processing unit and control unit of each element must be increased, which increases manufacturing cost. For these reasons many practical measurement apparatuses based on the PAT technique detect an acoustic wave from an object in a specific direction using a limited sized probe.
An example of such an apparatus is the PAT of a plate type measurement system disclosed in PTL 1. According to PTL 1, light is irradiated onto an object sandwiched by plates, and an acoustic wave is detected by an acoustic wave detector installed on the plate. In some cases the light is irradiated and the acoustic wave is detected for a plurality of times, and measured values are averaged, whereby such an effect as an improvement in the S/N ratio is implemented.
Image qualities (S/N ratio, artifact) of a reconstructed image in a photoacoustic diagnostic apparatus are influenced not only by the above mentioned acoustic wave detection conditions, but also by light irradiation conditions. If light is irradiated from outside an object, the level of light decays from the surface area to a deep part of the object due to the absorption of light by biological tissue. In effect it is difficult to irradiate light onto an object under ideal conditions.
It is possible to estimate the intensity of the irradiated light inside the object. In other words, the extent of decay of the level of irradiated light is determined considering an absorption coefficient according to a segment of the object, and light energy distribution is estimated. However it is difficult to completely eliminate the artifacts and errors of signal values even if the signal values of the photoacoustic wave are corrected based on the estimation result.
Furthermore, data becomes enormous if all factors that could influence the optical system in the imaging apparatus are recorded for the object. In the case of a method of recording and simulating the specifications and settings of each apparatus, the recording method and utilization method could be overly specific. It also involves complicated handling and time to calculate information on the irradiated light from the recorded data.
In the case of using a reconstructed image of the photoacoustic diagnostic apparatus for the purpose of medical diagnosis, it is necessary to know the degree of reliability and the influence of optical conditions within a reconstructed image. However conventionally photo-irradiation conditions depend on the design and specifications of the apparatus, and are recorded as apparatus-dependent information (e.g. setting values and imaging conditions of the apparatus). It is also time consuming to analyze information when conditions of the irradiated light are reproduced based on the recorded information. In photoacoustic diagnostic apparatuses as well, no technique is available to store information on individual irradiated light for each reconstructed image.
According to the technique disclosed in PTL 2, virtual light source information is stored to render computer graphics for creating a three-dimensional image. However information on irradiated light for imaging inside biological tissue cannot be stored as information corresponding to the reconstructed image. Therefore the only way to display and analyze a reconstructed image considering the conditions of irradiated light is to record redundant data which depends on the apparatus, and to perform time consuming analysis processing.